#include <kuai/misc/Matrix.hpp>
#include <kuai/mol/molalgo.hpp>
#include <kuai/mol/sss.hpp>

namespace kuai {
	bool sss(PhaseSpace& spaceQ, PhaseSpace& spaceP) {
		ConnectTable& ctQ = static_cast<PropertyConnectTable&>(*(spaceQ.find(MOL_INFO_CT)->second)).x_ct();
		const Index* prankQ = NULL;
		if (spaceQ.find(MOL_INFO_BOND_RANK) == spaceQ.end()) {
			ArrayProperty<Index>* p = new ArrayProperty<Index>(ctQ.countBonds());
			ConnectTable::RingFinder::getRingRankOfBonds(ctQ, p->x_ptr(), Index(-1));
			spaceQ.insert(MOL_INFO_BOND_RANK, SimObjPtr(p));
			prankQ = p->x_ptr();
		}
		else {
			prankQ = static_cast<PropertyIndex&>(*spaceQ[MOL_INFO_BOND_RANK]).x_ptr();
		}

		ConnectTable& ctP = static_cast<PropertyConnectTable&>(*(spaceP.find(MOL_INFO_CT)->second)).x_ct();
		const Index* prankP = NULL;
		if (spaceP.find(MOL_INFO_BOND_RANK) == spaceP.end()) {
			ArrayProperty<Index>* p = new ArrayProperty<Index>(ctP.countBonds());
			ConnectTable::RingFinder::getRingRankOfBonds(ctP, p->x_ptr(), Index(-1));
			spaceP.insert(MOL_INFO_BOND_RANK, SimObjPtr(p));
			prankP = p->x_ptr();
		}
		else {
			prankP = static_cast<PropertyIndex&>(*spaceP[MOL_INFO_BOND_RANK]).x_ptr();
		}

		return sss(ctQ, prankQ, ctP, prankP);
	}

	/*
	namespace { 
		class SubStrcutureSearch { 
		public: 
			SubStrcutureSearch(
				ConnectTable& query, const Index brQuery[],
				ConnectTable& pattern,  const Index brPattern[])
				: _q(&query), _p(&pattern)
			{
				nQA = _q->countAtoms();
				nPA = _p->countAtoms();
				assert(nQA >= nPA);

				nQB = _q->countBonds();
				nPB = _p->countBonds();
				assert (nQB >= nPB);

				_makeMappingMatrix(query, brQuery, pattern, brPattern, _aamap, _bbmap);

				_succeed = _makeMapping();
			}

		public:
			bool succeed() const { 
				return _succeed;
			};


		private:
			static void _makeMappingMatrix(
				ConnectTable& query, const Index brQuery[],
				ConnectTable& pattern,  const Index brPattern[],
				Matrix<char>& atomMap, Matrix<char>& bondMap) 
			{ 
				Index nQuery = query.countAtoms();
				Index nPattern = pattern.countAtoms();

				atomMap.reset(nQuery, nPattern, char(0));
				for (Index i = 0; i < nQuery; ++i) { 
					const ConnectTable::AtomType& atomI = query.getAtom(i);
					for (Index j = 0; j < nPattern; ++j) { 
						const ConnectTable::AtomType& atomJ = pattern.getAtom(j);
						if (atomI == atomJ && query.degree(atomI) >= pattern.degree(atomJ)) { 
							atomMap[i][j] = char(1);
						}
					}
				}

				nQuery = query.countBonds();
				nPattern = pattern.countBonds();
				bondMap.reset(nQuery, nPattern, char(0));
				for (Index i = 0; i < nQuery; ++i) { 
					const ConnectTable::BondType& bondI = query.getBond(i);
					int aI1 = query.indexOf(query.getAtom(bondI, 0));
					int aI2 = query.indexOf(query.getAtom(bondI, 1));
					for (Index j = 0; j < nPattern; ++j) { 
						const ConnectTable::BondType& bondJ = pattern.getBond(j);
						if (bondI != bondJ) 
						{
							continue;
						}

						int aJ1 = pattern.indexOf(pattern.getAtom(bondJ, 0));
						int aJ2 = pattern.indexOf(pattern.getAtom(bondJ, 1));

						if (brPattern[j] == 0 || brPattern[j] >= brQuery[i]) {
							if (atomMap[aI1][aJ1] && atomMap[aI2][aJ2] || atomMap[aI1][aJ2] && atomMap[aI2][aJ1]) { 
								bondMap[i][j] = char(1);
							}
						}
					}
				}
			}

			bool _makeMapping() { 
				std::vector<Index> countPossible(nPA, 0);
				for (Index i = 0; i < nQA; ++i) { 
					for (Index j = 0; j < nPA; ++j) { 
						if (aamap(i, j)) { 
							countPossible[j] += 1;
						}
					}
				}
				std::vector<Index>::const_iterator it = std::min_element(countPossible.begin(), countPossible.end());
				if (*it > 0) { 
					_aaIndex = _bbIndex = 0;
					_atomIndexQ.resize(_q->countAtoms(), 0);
					_atomIndexP.resize(_p->countAtoms(), 0);
					_bondIndexQ.resize(_q->countBonds(), 0);
					_bondIndexP.resize(_p->countBonds(), 0);
					Index startP = it-countPossible.begin();
					for (Index i = 0; i < nQA; ++i) { 
						if (aamap(i, startP)) { 
							if (_travel(i, startP, -1, -1)) { 
								return true;
							}
						}
					}
				}
				return false;
			}

			bool _travel(Index startQ, Index startP, Index fromQ, Index fromP) { 
				assert (aamap(startQ, startP));
				Index oldAAIndex = _aaIndex, oldBBIndex = _bbIndex;

				_mapAtom(startQ, startP);

				ConnectTable::AtomType& atomP = _p->getAtom(startP);
				ConnectTable::AtomType& atomQ = _q->getAtom(startQ);

				Index degreeP = _p->degree(atomP);
				Index degreeQ = _q->degree(atomQ);

				std::vector<Index> indexes(degreeP, 0);
				for (Index i = 0; i < degreeP; ++i) { 
					ConnectTable::AtomType& atomPI = _p->getNeighborAtom(atomP, i);
					Index iap = _p->indexOf(atomPI);
					if (iap == fromP) { 
						continue;
					}
					ConnectTable::BondType& bondPI = _p->getNeighborBond(atomP, i);
					Index ibp = _p->indexOf(bondPI);
					if (_bondIndexP[ibp] == 0) { 
						Index j = 0;
						for (j = indexes[i]; j < degreeQ; ++j) { 
							ConnectTable::AtomType& atomQJ = _q->getNeighborAtom(atomQ, j);
							ConnectTable::BondType& bondQJ = _q->getNeighborBond(atomQ, j);
							Index iaq = _q->indexOf(atomQJ);
							Index ibq = _q->indexOf(bondQJ);

							if (_atomIndexP[iap] == _atomIndexQ[iaq]
							&& _bondIndexP[ibp] == _bondIndexQ[ibq]
							&& aamap(iaq, iap) && bbmap(ibq, ibp))
							{
								_mapBond(bondQJ, bondPI);
								if (_atomIndexP[iap] > 0) 
								{ 
									break;
								}
								else if (_travel(iaq, iap, startQ, startP)) {
									break;
								}
								else { 
									_restoreBond();
									continue;
								}
							}
						}	// for (j = 0; j < degreeQ; ++j)

						if (j == degreeQ) { 
							_restore(oldAAIndex, oldBBIndex);
							Index k = i-1;
							for (; k < degreeP; --k) { 
								++indexes[k];
								if (indexes[k] < degreeQ) { 
									break;
								}
							}
							if (k == Index(-1)) { 
								return false;
							}
							for (k = k+1; 0 < k && k < degreeP; ++k) { 
								indexes[k] = 0;
							}
							i=-1;

							_mapAtom(startQ, startP);
						}
					}	// if (_bondIndexP[_p->indexOf(bondPI)]  == 0)
				}	// for (Index i = 0; i < degreeP; ++i)

				return true;
			}

			Index _mapAtom(Index iQ, Index iP) { 
				assert (_atomIndexQ[iQ] == 0);
				assert (_atomIndexP[iP] == 0);

				_atomIndexQ[iQ] = _atomIndexP[iP] = ++_aaIndex;
				return _aaIndex-1;
			}

			Index _mapBond(ConnectTable::BondType& bondQ, ConnectTable::BondType& bondP) { 
				Index iQ = _q->indexOf(bondQ), iP = _p->indexOf(bondP);
				assert (_bondIndexQ[iQ] == 0);
				assert (_bondIndexP[iP] == 0);

				_bondIndexQ[iQ] = _bondIndexP[iP] = ++_bbIndex;
				return _bbIndex-1;
			}

			void _restore(Index oldAAIndex, Index oldBBIndex) { 
				_restore(_atomIndexQ, oldAAIndex);
				_restore(_atomIndexP, oldAAIndex);
				_restore(_bondIndexQ, oldBBIndex);
				_restore(_bondIndexP, oldBBIndex);
				_aaIndex = oldAAIndex;
				_bbIndex = oldBBIndex;
			}

			void _restore(std::vector<Index>& array, Index oldIndex) { 
				for (std::vector<Index>::iterator 
					i = array.begin(); i != array.end(); ++i)
				{
					if (*i > oldIndex) { 
						*i = 0;
					}
				}
			}

			void _restoreBond() { 
				_bbIndex = _bbIndex - 1;
				_restore(_bondIndexQ, _bbIndex);
				_restore(_bondIndexP, _bbIndex);
			}

		public:
			SubMolecule<ConnectTable> getResult() { 
				assert (_aaIndex >= _p->countAtoms() && _bbIndex >= _p->countBonds());
				std::vector<ConnectTable::AtomType*> i2a(_aaIndex+1, NULL);
				for (Index i = 0; i < _atomIndexQ.size(); ++i)
				{ 
					if (_atomIndexQ[i] != 0) { 
						assert (_atomIndexQ[i] < i2a.size());
						assert (i2a[_atomIndexQ[i]] == NULL);
						i2a[_atomIndexQ[i]] = &(_q->getAtom(i));
					}
				}
				std::vector<ConnectTable::AtomType*> atoms; atoms.reserve(_p->countAtoms());
				std::vector<ConnectTable::BondType*> bonds; bonds.reserve(_p->countBonds());

				for (Index i = 0; i < _atomIndexP.size(); ++i) { 
					atoms.push_back(i2a[_atomIndexP[i]]);
				}
				Index nBonds = _p->countBonds();
				for (Index i = 0; i < nBonds; ++i) { 
					ConnectTable::BondType& bondI = _p->getBond(i);
					ConnectTable::AtomType& atom1 = _p->getAtom(bondI, 0);
					ConnectTable::AtomType& atom2 = _p->getAtom(bondI, 1);
					Index index1 = _p->indexOf(atom1);
					Index index2 = _p->indexOf(atom2);
					index1 = _atomIndexP[index1];
					index2 = _atomIndexP[index2];
					bonds.push_back(_q->getBond(*i2a[index1], *i2a[index2]));
				}

				return SubMolecule<ConnectTable>(*_q, atoms, bonds);                
			}

		private:
			ConnectTable* _q;
			ConnectTable* _p;
			Matrix<char> _aamap;
			Matrix<char> _bbmap;

			bool aamap(Index i, Index j) { 
				return _aamap[i][j] != 0;
			}
			bool bbmap(Index i, Index j) { 
				return _bbmap[i][j] != 0;
			}

			Index _aaIndex, _bbIndex;
			std::vector<Index> _atomIndexQ, _atomIndexP, _bondIndexQ, _bondIndexP;
			Index nQA, nPA, nQB, nPB;
			bool _succeed;
		};
	}


	bool sss(ConnectTable& ctQ, const Index* bondRankQ, ConnectTable& ctP, const Index* bondRankP) {
		if (ctQ.countAtoms() < ctP.countAtoms() || ctQ.countBonds() < ctP.countBonds()) {
			return false;
		}
		std::vector<Index> vbondRankQ;
		if (bondRankQ == NULL) {
			vbondRankQ.resize(ctQ.countBonds());
			RingFinder<ConnectTable>::getRingRankOfBonds(ctQ, &vbondRankQ[0], ctP.countBonds());
			bondRankQ = &vbondRankQ[0];
		}

		std::vector<Index> vbondRankP;
		if (bondRankP == NULL) {
			vbondRankP.resize(ctP.countBonds());
			RingFinder<ConnectTable>::getRingRankOfBonds(ctP, &vbondRankP[0], ctP.countBonds());
			bondRankP = &vbondRankP[0];
		}

		SubStrcutureSearch search(ctQ, bondRankQ, ctP, bondRankP);
		return search.succeed();
	}

	/*

	namespace {
	const Index MAX_RING_RANK = 6;

	inline Char* concat(const Char* source, Char* dest) {
	while (*source) {
	*dest = *source;
	++dest, ++source;
	}
	return dest;
	}

	}

	SubStructureSearchPattern::SubStructureSearchPattern(const Index data[]) 
	: ConnectTable(data)
	{
	_ringrank = data + ConnectTable::count_size(count_atoms(), count_bonds());
	}

	SubStructureSearchPattern::SubStructureSearchPattern(Index nAtoms, const Index atoms[], Index nBonds, const Index bonds[], Index buffer[])  
	: ConnectTable(nAtoms, atoms, nBonds, bonds, buffer)
	{
	Index* pring = buffer + ConnectTable::count_size(nAtoms, nBonds);
	std::fill(pring, pring+nBonds, INVALID_INDEX);
	for (Index i = 0; i < nBonds; ++i) {
	if (pring[i] == INVALID_INDEX) {
	_setup_ringrank(i, MAX_RING_RANK);
	}
	else if (pring[i] > 3){
	_setup_ringrank(i, pring[i]-1);
	}
	}

	_ringrank = pring;
	}

	Index SubStructureSearchPattern::count_size(Index nAtoms, Index nBonds) {
	return ConnectTable::count_size(nAtoms, nBonds) + nBonds;
	}

	void SubStructureSearchPattern::_setup_ringrank(Index ibond, Index maxring) const {
	// TODO
	}

	void SubStructureSearchPattern::setup_filter(Index version, Byte result[], Index n) const 
	{
	Char buffer[64];
	for (Index i = 0; i < count_atoms(); ++i) {
	Index v = _filter_atom(i, buffer) % n;
	result[v/8] |= BIT_FLAGS[v%8];

	Index d = degree(i);
	for (Index j = 1; j < d; ++j) {
	Index a1 = neighbor_atom(i, j);
	Index bi1 = neighbor_bond(i, j);
	for (Index k = 0; k < j; ++k) {
	Index a2 = neighbor_atom(i, k);
	Index bi2 = neighbor_bond(i, k);
	v = _filter_angle(a1, bi1, i, bi2, a2, buffer) % n;
	result[v/8] |= BIT_FLAGS[v%8];
	}
	}
	if (d == 3) {
	Index a1 = neighbor_atom(i, 0), b1 = neighbor_bond(i, 0);
	Index a2 = neighbor_atom(i, 1), b2 = neighbor_bond(i, 1);
	Index a3 = neighbor_atom(i, 2), b3 = neighbor_bond(i, 2);
	v = _filter_fork(i, a1, b1, a2, b2, a3, b3, buffer) % n;
	result[v/8] |= BIT_FLAGS[v%8];
	}
	}
	for (Index i = 0; i < count_bonds(); ++i) {
	Index v = _filter_bond(i, buffer) % n;
	result[v/8] |= BIT_FLAGS[v%8];

	Index a1 = bond_atom_1(i);
	Index a2 = bond_atom_2(i);
	Index d1 = degree(a1);
	Index d2 = degree(a2);

	for (Index j = 0; j < d1; ++j) {
	Index a0 = neighbor_atom(a1, j);
	if (a0 != a2) {
	Index b01 = neighbor_bond(a1, j);
	for (Index k = 0; k < d2; ++k) {
	Index a3 = neighbor_atom(a2, k);
	if (a3 != a1 && a3 != a0) {
	Index b23 = neighbor_bond(a2, k);
	v = _filter_torsion(a0, b01, a1, i, a2, b23, a3, buffer) % n;
	result[v/8] |= BIT_FLAGS[v%8];
	}
	}
	}
	}
	}
	}

	Index SubStructureSearchPattern::_hash(const Char* p) 
	{
	size_t result = 2166136261U;
	for(; *p; ++p) {
	result = 16777619U * result ^ size_t(*p);
	}
	return Index(result);
	}

	Index SubStructureSearchPattern::_filter_atom(Index i, Char buffer[]) const {
	return _hash(symbol(i).c_str());
	}

	Index SubStructureSearchPattern::_filter_bond(Index i, Char buffer[]) const {
	Index a1 = bond_atom_1(i);
	Index a2 = bond_atom_2(i);
	if (atom(a1) > atom(a2)) {
	std::swap(a1, a2);
	}

	Char* p = buffer;
	p = concat(symbol(a1).c_str(), p);
	*p++ = BOND_SYMBOLS[bond_order(i)];
	p = concat(symbol(a2).c_str(), p);
	*p = NULL;
	return _hash(p);
	}

	Index SubStructureSearchPattern::_filter_angle(Index i, Index bij, Index j, Index bjk, Index k, Char buffer[]) const {
	if (atom(i) > atom(j)) {
	std::swap(i, j);
	std::swap(bij, bjk);
	}

	Char* p = buffer;
	p = concat(symbol(i).c_str(), p);
	*p++ = BOND_SYMBOLS[bond_order(bij)];
	p = concat(symbol(j).c_str(), p);
	*p++ = BOND_SYMBOLS[bond_order(bjk)];
	p = concat(symbol(k).c_str(), p);
	*p = NULL;
	return _hash(p);
	}
	Index SubStructureSearchPattern::_filter_torsion(Index i, Index bij, Index j, Index bjk, Index k, Index bkl, Index l, Char buffer[]) const {
	if (atom(j) == atom(k) && atom(i) > atom(l) || atom(j) > atom(k)) {
	std::swap(i, l); 
	std::swap(j, k);
	std::swap(bij, bkl);
	}
	Char* p = buffer;
	p = concat(symbol(i).c_str(), p);
	*p++ = BOND_SYMBOLS[bond_order(bij)];
	p = concat(symbol(j).c_str(), p);
	*p++ = BOND_SYMBOLS[bond_order(bjk)];
	p = concat(symbol(k).c_str(), p);
	*p++ = BOND_SYMBOLS[bond_order(bkl)];
	p = concat(symbol(l).c_str(), p);
	*p = NULL;
	return _hash(p);
	}
	Index SubStructureSearchPattern::_filter_fork(Index center, Index a1, Index b1, Index a2, Index b2, Index a3, Index b3, Char buffer[]) const {
	if (atom(a1) > atom(a2)) {
	std::swap(a1, a2);
	std::swap(b1, b2);
	}
	if (atom(a2) > atom(a3)) {
	std::swap(a2, a3);
	std::swap(b2, b3);
	}
	if (atom(a1) > atom(a2)) {
	std::swap(a1, a2);
	std::swap(b1, b2);
	}

	Char* p = buffer;
	p = concat(symbol(center).c_str(), p);
	*p++ = '/';
	*p++ = BOND_SYMBOLS[bond_order(b1)];
	p = concat(symbol(a1).c_str(), p);
	*p++ = ',';
	*p++ = BOND_SYMBOLS[bond_order(b2)];
	p = concat(symbol(a2).c_str(), p);
	*p++ = ',';
	*p++ = BOND_SYMBOLS[bond_order(b3)];
	p = concat(symbol(a3).c_str(), p);
	*p = NULL;
	return _hash(p);
	}

	*/

}
